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ABSTRACT 

Recently,  a  method  for  generating  an  ARMA 
spectral  estimator  model  which  possessed  super- 
resolution  performance  was  developed  [2].  This 
method  entailed  minimizing  a  weighted  quadratic 
functional  of  a  set  of  "basic  error  terms.’*  An 
issue  which  remained  to  be  resolved  at  that  time 
was  the  selection  of  the  weighting  matrix  that 
characterized  the  functional  being  minimized.  A 
weighting  matrix  selection  procedure  has  recently 
been  developed  and  is  herein  reported  [  8 ] .  This 
orocedure  has  typically  yielded  an  improvement  in 
spectral  estimation  performance. 

I.  INTRODUCTION 

In  this  paper  we  shall  be  concerned  with  the 
cask  of  estimating  the  power  spectral  density  of  a 
zero  mean,  wide  sense  stationary  random  time  series 
{x(n)}  from  a  finite  set  of  observations.  To  this 
end,  knowledge  of  the  time  series*  underlying  auto* 
correlation  sequence  as  formally  defined  by 

rx(n)  *  E(x(n+k)  x*(k);  (1) 

conveys  all  the  information  required.  Here,  E(  } 
denotes  the  expected  value  operators  and  *  denotes 
the  operation  of  complex  conjugation.  The  time 
series  is  characterized  in  the  frequency  domain  by 
its  power  spectral  density  as  given  by 

S8u>  *  7  rx(n)e_j"n  (2) 

n»— » 

which  is  recognized  as  being  the  Fourier  transform 
of  the  autocorrelation  sequence. 

Upon  examination  of  (1)  and  (2)  it  is  apparent 
chat  determination  of  the  time  series  power  spectral 
density  entails  complete  knowledge  of  the 
generally  infinite  length  autocorrelation  sequence. 
Here,  we  will  be  concerned  with  extracting  this 
information  from  the  finite  sec  of  time  series 
observations 

x(i).  x(2) . x(N)  *  (3) 
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Unless  some  constraints  are  imposed  on  the 
time  series'  basic  nature,  however,  there  exists 
a  fundamental  incompatibility  in  estimating  the 
required  statistical  knowledge  from  the  finite  sec 
of  data.  This  dilemma  is  usually  resolved  by 
postulating  a  finite  parameter  linear  model  to 
represent  the  time  series.  In  terms  of  parameter 
parsimony,  the  causal  autoregressive  moving  average 
(ARMA)  model  of  order  (p,q)  as  specified  by 


+  l  a.x(n-i)  *  l  b.e(n-j) 
i-1  i*0  J 


is  generally  the  most  effective  linear  model  [1], 

In  this  model,  the  (unobserved)  excitation  process 
(e(n)}  is  assumed  to  be  zero-mean,  unit  variance 
Gaussian  whits  noise.  It  is  important  to  note  that 
the  more  specialized  autoregressive  model 
(i.e.,  bj  =  0  for  J  4  0)  generally  requires  a  much 
higher  model  order  p  to  achieve  comparable 
spectral  estlmacss.  Conceptually,  then,  the  more 
general  ARMA  model  is  the  logical  model  choice. 

It  is  well  known  that  the  power  spectral  dan- 
slty  of  a  process (x(n)} chat  satisfies  (4)  is  given 

*  11  +  a.e-^...  +a  rjp“|2  <3> 

I  1  P  f 

Thus  the  cask  of  estimating  the  power  spectral 
density  of  the  time  series  can  be  accomplished  by 
estimating  the  ARMA  model  parameters  a*  and  b j . 

Several  procedures  for  estimating  the  a*  and 
bj  parameters  have  recently  been  developed  [2-8]. 

Of  these  procedures,  one  developed  by  Cadzow  [2] 
has  been  shown  to  be  effective  in  a  variety  of 
cases.  The  crux  of  this  procedure  lies  In  obtain¬ 
ing  the  autoregressive  parameters  by  minimizing  a 
weighted  quadratic  function  of  a  set  of  zero  mean 
error  elements.  It  was  pointed  out  in  [2]  chat  the 
ef feccivenesa  of  this  procedure  is  dependent  on  a 
judicious  selection  of  the  weighting  elements  in 
the  quadratic  function.  This  papar  develops  an 
alternative  weighting  element  selection  to  chat 
used  in  [2]  which  results  in  improved  spectral 
estimates.  _  t  ■  i  " 
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II.  THE  MODEL  EQUATION  ERROR 
SPECTRAL  ESTIMATOR. 


and 


The  spectral  estimation  procedure  of  this 
paper  is  predicated  on  the  procedure  in  [2].  For 
completeness,  this  procedure  Is  discussed  below. 


Of  primary  Importance  in  spectral  estimation 
Is  the  method  for  estimating  the  autoregressive  co¬ 
efficients  a^  in  equation  (4).  An  effective  method 
for  estimating  these  coefficients  entails  multiply¬ 
ing  both  sides  of  (4)  by  the  term  x*(n-m)  to  yield 
the  "basic  error  terms” 


e(m,n) 


P 

-  x(n)x*(n-m)  +  £  a1x(n-i)x*(n-m)  (6a) 

1-1 


*  7  b (n-j )x*(n-m) 

j-0  2 


for  q+1  ^m  <_N-1 
max  ( p+1 ,  m+1 )  <n  <^N 


(6b) 


where  the  range  on  the  m  and  n  variables  Is  dictat¬ 
ed  by  the  time  series  observation  range  1  <^N. 

If  the  time  series  Is  in  fact  an  ARMA  process  of 
order  less  than  or  equal  to  (p,q)  Chen  the  basic 
error  terms  are  each  zero  mean  random  variables. 
Furthermore,  the  basic  error  terms  are  seen  to  be 
functions  of  the  autoregressive  coefficients 
al»  a2»  ...»  ap.  With  these  two  properties  in  mind, 
a  reasonable  selection  of  the  autoregressive  co¬ 
efficients  is  one  chat  causes  each  of  the  e(m,n) 
terms  to  be  as  close  to  its  mean  value  of  zero  as 
possible.  This  goal  is  achieved  by  minimizing  a 
squared-error  criterion  of  the  form 


f(a)  * 


(7) 


where  e  Is  a  vector  of  appropriately  arranged 
e(m,n)  terms,  W  is  a  positive  semidefinite  weight¬ 
ing  matrix,  and  *  denotes  the  operation  of  complex 
conjugate  transposition. 


A  more  specific  format  of  the  minimization 
criterion  (7),  and  the  one  considered  in  [21,  is 


N*-i  i  N  2 

:  (a)  -  ;  w(m) ;  2  •(«,«)■ 

m»q+l  n-s 


s  -  max(p+l ,uH“l ) 

(8) 


N-i  N  N  * 

c..  •  III  v(m)x(n-k)xU-m)x  (n-m)x  (2-i) 
rn-q+1  h-s  I- s 

s  *  maximal ,  p+1 )  ( lQd) 

An  improvemant  in  the  above  autoregressive  co¬ 
efficient  estimation  procedure  can  be  realized  by 
also  considering  the  backward  varsion  of  the  time 
series  (x(n)h  Also,  an  estimate  of  the  numerator 
spectrum  B(w)  in  (5)  must  be  obtained  in  order  to 
arrive  at  the  complete  power  spectral  density  esti¬ 
mate  of  (x(n)}.  The  details  of  these  two  tasks  are 
presented  in  [2]. 


III.  WEIGHTING  ELEMENT  SELECTION 


In  order  to  obtain  autoregressive  parameters 
using  the  above  procedure,  the  elements  v(m)  in  (8) 
must  first  be  selected.  In  [21  the  weights 

w(m)  -  (N-m)“*  q+1  f.  m  ^  M-l  (11) 

were  employed. 


A  more  prudent  weight  selection  can  be 
developed  by  considering  the  random  terms 


Nl  2 

l  e(m,n) 

n-s 


q+1  i  a  1  N-l 


(12) 


associated  with  the  weights  v(m)  in  (8)  [ 3  ] .  A 
logical  selaction  for  the  w(m)  weights  is  the 
inverse  of  the  variances  of  the  terms  in  (12), 
that  is 


w(m)  •  \ 

var 

“  N  ,2] 

l  e(o,n) |  i 

i 

J  n-s  1  J 

(13) 


In  this  way,  the  terms  in  the  minimization  criter¬ 
ion  (8)  which  have  smaller  variances  from  their 
mean  valua  are  weighted  proportionately  higher  than 
those  terms  with  larger  variances  from  their  mean. 


In  chis  expression  the  w(m)  are  nonnegative 
weighting  elements. 


Using  standard  calculus,  it  is  readily  shown 
that  the  sec  of  autoregressive  coefficients  which 
minimize  (8)  are  given  by: 


a  - 

c-1c 

(9) 

where 

a  - 

[ a  ] 

(10a) 

C  • 

uik  ?'? 

(10b) 

^  m 

-(cIQ  c:o  ...  cp0; 

(10c) 

It  is  easily  shown 

chat 

N 

!  21 

N 

var 

I  7  e(m,a) 

j-  - 

7  r  ( l-n)  cU-n) 

1 

Jn-s 

'  J  n-s  j 

>s 

s  ■  max(art*l,p+i) 

wh,r#  [  a-m 

1  l  b 

!  i-0 

ibL»  , 

0  2  ®  i  q 

C(B)  c*(  ini )  . 

1 

-1  _  m  -q 

i 

»  0 
i 

• 

otherwise 

(14a) 


(14b) 


Unfortunately,  tha  desired  variances  are  seen 
to  be  dependant  on  the  unknown  parameters 
bo*  b^>  . ...  bq.  However,  an  approximate  express¬ 
ion  for  the  inverse  variance  weights  of  (13)  can 
be  realized  if  a  reasonable  approximation  of  the 
c(m)  elements  in  (14b)  can  be  found. 

One  can  gain  Insight  about  the  structure  of 
the  c(a)  elements  by  forming  the  polynomials  B(z) 
and  C(z)  defined  as 

B(z)  -  brtzq  +  b.  z^"1  +  . . .  +b  .z+b  (15a) 

0  l  q-i  q 

C(z)  •  c^q)!**1  +  c(-q+l)z'<T+i  +  . ..  +c(0)  +c(l)z 

+  ...  +c(q)zq 

(15b) 

It  is  easily  shown  chat 

C(z)  •  B(z)B(z*1)  (16) 


3t  -  8^  i  -  2,4,6,  ....  q  (20b) 

where  ImtS^]  denotes  the  imaginary  part  of  3*. 

If  q  is  odd  the  unpaired  zero  is  assumed  to  be 
uniformly  distributed  on  the  real  axis  Inside  the 
unit  circle, 

t  ,  *1  <  J,  *  l  ttUJ  *  0  <20c) 

fg  (8q)  -|2  “  q“  ’ 

H  0  ,  otherwise 

Using  these  assumptions  about  the  zeroes  of 
B(z),  one  can  straightforwardly  calculate  the 
desired  approximation  to  C(z)  by  determining  the 
expected  value  of  expression  (18) 

d 

C(*)  -  EU"*1  B  (z-a.)d-ze,)}  (21) 

i-0 

By  carrying  out  this  calculation,  one  finds  that 
for  a  complex  time  series  (x(n)} 


Furthermore,  B(z)  can  be  factored  as 


B(z )  -  bQ  if  (z-it) 
u  i-i 


where  the  3*  are  zeroes  of  the  polynomial  B(z). 
Applying  (16)  it  is  found  that 

C(z)  «bn2z'q  3  (z-3.)(l-zSt)  (18) 

0  1-1  1 

2 

Thus,  C(z)  can  be  found  (to  within  the  constant  b(f) 
using  (18)  from  knowledge  of  Che  q  zeroes 
3^  or  B(z) . 

A  reasonable  approximation  of  C(z)  and  there¬ 
fore  of  the  c(m)  elements  can  be  found  by 
approximating  the  location  of  the  zeroes  of  3(z). 
One  3uch  approximation  is  realized  by  assuming  chat 
each  zero  is  a  random  variable  uniformly  distribut¬ 
ed  within  the  complex  unit  circle,1-  so  that  its 
probability  density  function  is 


Mil 

i  ■  1,2, . . .  ,q 

otherwise 


(19) 


If  the  time  series  (x(n);  is  a  real  process, 
then  the  zeroes  of  B(z)  must  form  complex  conjugate 
pairs.  For  this  case  it  is  assumed  chat  3/2  of  the 
zeroes  are  uniformly  distributed  within  the  upper 
half  of  the  complex  unit  circle,  chat  is 


Cc(z)  -  1 

and  for  a  real  time  series 

L  "J 

where 

i  ,  q  even 

k  - 

^2~  .  q  odd 


(22a) 

(22b) 


Thus,  the  approximate  inverse  variance  weights  are 
given  by 

r  n  n  ^  v1 

w(m)  -  7  [  r  (l-n)c(l-n)  I  q+1  <  m  <  N-L 


where  s  -  max(nH-l,p+l)  and  the  c(m)  elements  are 
the  coefficients  corresponding  to  the  za  terms  of 
the  polynomials  (22a)  or  (22b). 


IV.  NUMERICAL  EXAMPLE 

In  order  to  compare  the  effectiveness  of  the 
new  ARMA  spectral  estimator  with  the  estimator 
In  [2],  Che  classical  problem  of  resolving 
two  closely  spaced  (in  frequency)  sinusoids  in 
white  noise  will  be  considered.  Specifically,  the 
time  series  under  study  is  specified  by 


it;  ^  l  and  ImUjJ  ^  0  (20a) 

otherwise 

i  -  1.3, - q-i 


a  zero  of  B(z)  is  outside  the  unit  circle,  then 
cue  corresponding  zero  of  Su*1)  is  inside  the  unit 
circle,  and  from  equation  (16)  it  is  clear  that 
c;z)  will  not  be  affected. 


x(n)  •  f'lo  cos(0.4irn)  *  *2  cos(0.426wn)  +w(n)  (24) 

where  Cw(n);  is  a  white  Gaussian  noise  process  of 
zero  mean  and  unit  variance.  The  sinusoids  of 
normalized  frequencies  0.4  and  0.426  are  readily 
calculated  to  have  signal -co-noise  ratios  (SNR)  of 
lOdB  and  OdB,  respectively.  A  sequence  of  length 
512  defined  over  0  <  n  <  511  was  next  generated 
using  this  relationship.  Furthermore,  in  order  to 
provide  a  statistical  basis  for  our  comparison. 


this  512  length  sequence  wee  then  decomposed  Into 
eight  disjoint  sequences  each  of  length  64  defined 
on  0  _<  n  63,  64  n  127,  . . . ,  448  _<  n  <_  511. 

An  ensemble  consisting  of  eight  subsequences  each 
of  length  64  has  thereby  been  generated  with  each 
subsequence  having  a  different  noise  sample  and  a 
different  initial  phase  between  the  two  sinusoids. 
This  latter  condition  is  useful  in  revealing  any 
potential  sensitivity  to  initial  phase  chat  the 
new  ARMA  spectral  estimation  method  may  possess. 

The  spectral  estimates  which  resulted  when 
the  (N-m)^  weights  and  the  new  Inverse  variance 
weights  were  applied  to  the  ARMA  spectral  estima¬ 
tor  are  displayed  in  Figures  la  and  lb,  respective* 
ly.  The  ordinates  are  scaled  from  -20dB  to  60dB  for 
each  individual  plot.  In  both  cases  the  spectral 
estimator  ordar  was  (15,15).  It  is  clear  that  the 
inverse  variance  weight  estimator  was  abla  to 
resolve  the  two  sinusoids  in  more  cases  than  the 
(N-m)^  weight  estimator  could.  Moreover,  the 
incidence  of  false  peaks  in  the  inverse  variance 
weight  estimates  la  smaller  than  that  of  the 
(N-m)4  weight  estimates. 


V.  CONCLUSIONS 


An  improved  weight  selection  for  a  recently 
*  developed  ARMA  spectral  estimation  procedure  was 

developed  .^The  autoregressive  parameters  are 
found  in  this  procedure  by  minimizing  a  weighted 
\  sum  of  squares  of  zero  mean  basic  error  terms. 

i  The  new  weight  selection  is  chosen  to  provide  more 

;  heavy  weighting  to  those  terms  la  the  sum  which 

A  possess  lower  variances.  Empirical  evidence 

t  indicates  chat  this  new  weight  selection  provides 

superior  spectral  estimation  performancs  when  com¬ 
pared  to  the  original  weight  selection. 
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Figure  1.  Spectral  Estimates  of  Two  Sinusoids  at 
Normalized  Frequencies  of  O.iOO(lOdB)  and  0.426 
(OdB)  in  Additive  White  Noise,  a)  (N-m)**  weights 
used,  b)  Inverse  variance  weights  used. 
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